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This paper describes the Queueing Network Analyzer (QNA), a software 
package developed at Bell Laboratories to calculate approximate congestion 
measures for a network of queues. The first version of QNA analyzes open 
networks of multiserver nodes with the first-come, first-served discipline and 
no capacity constraints. An important feature is that the external arrival 
processes need not be Poisson and the service-time distributions need not be 
exponential. Treating other kinds of variability is important. For example, 
with packet-switched communication networks we need to describe the conges- 
tion resulting from bursty traffic and the nearly constant service times of 
packets. The general approach in QNA is to approximately characterize the 
arrival processes by two or three parameters and then analyze the individual 
nodes separately. The first version of QNA uses two parameters to characterize 
the arrival processes and service times, one to describe the rate and the other 
to describe the variability. The nodes are then analyzed as standard GI/G/m 
queues partially characterized by the first two moments of the interarrival- 
time and service-time distributions. Congestion measures for the network as 
a whole are obtained by assuming as an approximation that the nodes are 
stochastically independent given the approximate flow parameters. 

I. INTRODUCTION AND SUMMARY 

Networks of queues have proven to be useful models to analyze the 
performance of complex systems such as computers, switching ma- 
chines, communications networks, and production job shops. 1-7 To 
facilitate the analysis of these models, several software packages have 
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been developed in recent years, e.g., BEST/1, 8 CADS, 9 PANA- 
CEA, 10-12 and one based on Heffes. 13 These software packages contain 
algorithms for Markov models that can be solved exactly. For some 
applications, the model assumptions are at least approximately satis- 
fied, so that the analysis can be very helpful. For many other appli- 
cations, however, the model assumptions are not even approximately 
satisfied, so that the analysis can be misleading. 

A natural alternative to an exact analysis of an approximate model 
is an approximate analysis of a more exact model. This paper describes 
a software package called the Queueing Network Analyzer (QNA), 
which was recently developed at Bell Laboratories to calculate ap- 
proximate congestion measures for networks of queues. QNA goes 
beyond existing exact methods by treating non-Markov networks: The 
arrival processes need not be Poisson and the service-time distribu- 
tions need not be exponential. QNA treats other kinds of variability 
by approximately characterizing each arrival process and each service- 
time distribution with a variability parameter. It is also possible to 
analyze large networks quickly with QNA because the required calcu- 
lations are minimal, the most complicated part being the solution of a 
system of linear equations. The current version of QNA is written in 
FORTRAN. 

Here is a rough description of the model: There is a network of 
nodes and directed arcs. The nodes represent service facilities and the 
arcs represent flows of customers, jobs, or packets. There is also one 
external node, which is not a service facility, representing the outside 
world. Customers enter the network on directed arcs from the external 
node to the internal nodes, move from node to node along the internal 
directed arcs, and eventually leave the system on one of the directed 
arcs from an internal node to the external node. The flows of customers 
on the arcs are assumed to be random so that they can be represented 
as stochastic processes. 

If all servers are busy at a node when a customer arrives, then the 
customer joins a queue and waits until a server is free. When there is 
a free server, that customer begins service, which is carried out without 
interruption. Successive service times at each node are assumed to be 
random variables, which may depend on the type of customer but 
which otherwise are independent of the history of the network and are 
mutually independent and identically distributed. After the customer 
completes service, he goes along some directed arc from that node to 
another node. The customer receives service in this way from several 
internal nodes and then eventually leaves the network. A picture of a 
typical network (without the external node) is given in Fig. 1. 

An important feature of the model is that there may be flows from 
node j to node i, as well as flows from node i to node j. This is of 
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Fig. 1 — An open network of queues. 



course useful when customers can return to a node where they previ- 
ously received service, but it is also useful when customers cannot 
return to a node where they previously received service. Then flows 
from node j to node i represent different customers than the customers 
that flow from node i to node j. 

To be precise about the model, we give a list of basic assumptions. 
It is worth noting, however, that work is under way to extend QNA so 
that it can analyze systems in which each of the following assumptions 
is replaced by obvious alternatives. The general approximation tech- 
nique is flexible, so that it is not difficult to modify and extend the 
algorithm. 

Assumption 1. The network is open rather than closed. Customers 
come from outside, receive service at one or more nodes, and eventually 
leave the system. 

Assumption 2. There are no capacity constraints. There is no limit 
on the number of customers that can be in the entire network and 
each service facility has unlimited waiting space. 

Assumption 3. There can be any number of servers at each node. 
They are identical independent servers, each serving one customer at 
a time. 

Assumption 4. Customers are selected for service at each facility 
according to the first-come, first-served discipline. 

Assumption 5. There can be any number of customer classes, but 
customers cannot change classes. Moreover, much of the analysis in 
QNA is done for the aggregate or typical customer (see Sections 2.3 
and VI). 

Assumption 6. Customers can be created or combined at the nodes, 
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e.g., an arrival can cause more than one departure (see Section 2.2). 
(Think of messages.) 

The general approach is to represent all the arrival processes and 
service-time distributions by a few parameters. The congestion at each 
facility is then described by approximate formulas that depend only 
on these parameters. The parameters for the internal flows are deter- 
mined by applying an elementary calculus that transforms the param- 
eters for each of the three basic network operations: superposition 
(merging), thinning (splitting), and flow through a queue (departure). 
These basic operations are depicted in Fig. 2. When the network is 
acyclic (e.g., queues in series), the basic transformations can be applied 
successively one at a time, but in general it is necessary to solve a 
system of equations or use an iterative method. To summarize, there 
are four key elements in this general approach: 

1. Parameters characterizing the flows and nodes that will be readily 
available in applications and that have considerable descriptive power 
in approximations of the congestion at each node. 

2. Approximations for multiseruer queues based on the partial infor- 
mation provided by the parameters characterizing the arrival process 
and the service-time distribution at each node. 

3. A calculus for transforming the parameters to represent the basic 
network operations: merging, splitting, and departure. 

4. A synthesis algorithm to solve the system of equations resulting 
from the basic calculus applied to the network. 

The current version of QNA uses two parameters to characterize 
the arrival processes and the service times, one to describe the rate 
and the other to describe the variability. (Three-parameter algorithms 
are being developed, however.) For the service times, the two param- 





(c) 

Fig. 2 — Basic network operations: (a) Superposition or merging, (b) Decomposition 
or splitting, (c) Departure or flow through a queue. 
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eters are the first two moments. However, we actually work with the 
mean service time r and the squared coefficient of variation c 2 , which 
is the variance of the service time divided by the square of its mean. 
The user has the option of working with the service rate n = t -1 
instead of t. For the arrival processes, the parameters are associated 
with renewal-process approximations. The first two parameters are 
equivalent to the first two moments of the renewal interval (interval 
between successive points) in the approximating renewal process. The 
equivalent parameters we use are the arrival rate X, which is the 
reciprocal of the renewal-interval mean, and the squared coefficient 
of variation c 2 , which is the variance of the renewal interval divided 
by the square of its mean. 

We obtain the approximation of the flows by applying the general 
framework and the basic procedures for approximating point processes 
in Whitt, 14 incorporating refinements such as the hybrid procedures 
developed for merging by Albin. 1516 Of course, the general idea of 
simple two-parameter approximations for stochastic point processes 
goes back at least to the equivalent random method for approximating 
overflow streams (see Wilkinson, 17 Cooper, 18 and references there). 
Renewal-process approximations for such point processes were intro- 
duced by Kuczura 19 (also see Rath and Sheng 20 ). Two-parameter 
approximations for networks of queues similar to QNA have also been 
developed by others, apparently first by Reiser and Kobayashi 21 (also 
see Kuehn, 22 Sevcik et al., 23 Chandy and Sauer, 24 Chapter 4 of Gelenbe 
and Mitrani, 4 and Shanthikumar and Buzacott 6 ). These two-parame- 
ter approximations for networks of queues are also similar in spirit to 
two -parameter approximations for networks of blocking systems with 
alternate routing (see Katz 25 ). 

Some authors have referred to these two-parameter heuristic ap- 
proximations for networks of queues as diffusion approximations, 4,21 
but diffusion processes are not actually used. Diffusion approximations 
and associated heavy-traffic limit theorems have motivated some of 
the heuristic approximations in the literature and in QNA, and they 
are closely related to the asymptotic method for approximating point 
processes, 14 but the heuristic approximations in QNA are not the same 
as the more complicated diffusion approximations for networks of 
queues in Iglehart and Whitt, 26 Harrison and Reiman, 27 and Rei- 
man. 28 ' 29 

The approximation method in QNA is perhaps best described as a 
parametric-decomposition method, 22 because the nodes are analyzed 
separately after the parameters for the internal flows are determined. 
Moreover, when the congestion measures are calculated for the net- 
work as a whole, the nodes are treated (approximately) as being 
stochastically independent. This independence can be interpreted as 
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a generalization of the product-form solution that is valid for Marko- 
vian networks, i.e., in the Markov models the components of the vector 
representing the equilibrium number of customers at each node are 
stochastically independent, so that the probability mass functions for 
the vector is the product of the probability mass functions for the 
components. While QNA can be thought of as a decomposition method 
or an extended-product-form solution, an effort is made to capture the 
dependence among the nodes. The idea is to represent this dependence 
approximately through the internal flow parameters. 

To see the motivation for QNA, consider the elementary open 
network containing a single node with a single server, an infinite 
waiting room, and the first-come, first-served discipline. Suppose there 
is a single customer class with each customer being served only once 
before departing. The standard Markov model of this elementary 
network, which is embodied in BEST/1, CADS, and PANACEA, is 
the classical M/M/l queue, 1,3,18 which has a Poisson arrival process 
and an exponential service-time distribution. For the M/M/l model, 
the expected waiting time EW (before the customer begins service) is 

EW = tp/(1 - p), (1) 

where r is the mean service time and p is the traffic intensity, which 
is assumed to satisfy * p < 1. 

On the other hand, QNA uses an approximation for the GI/G/1 
model to represent this one-node network. The GI/G/1 model has a 
renewal arrival process and both the interarrival-time distribution and 
the service-time distribution are general. In QNA, the arrival process 
is represented by a renewal process partially characterized by two 
parameters: the arrival rate X and the variability parameter c 2 . The 
service-time distribution is also partially characterized by two param- 
eters: the mean service time r and the variability parameter c 2 . In 
contrast to (1), the formula for the expected waiting time in QNA is 

EW = t P (cI + d)g/2(l - p), (2) 

where g = g(p, cl, cl) is either one (when c\ ^ 1) or less than one 
(when cl < 1); see (45). When g( Pt c 2 , cl) = 1, (2) differs from (1) by 
the factor (cl + c 2 )/2. When the arrival process is Poisson, cl = 1; 
when the service-time distribution is exponential, cl = 1. Hence, if the 
GI/G/1 model is actually an M/M/l model, (2) reduces to (1). Of 
course, the user of QNA can set cl = cl = 1 and obtain (1). In fact, 
the values cl = 1 and cl = 1 are default values that the program uses 
if the user does not have variability parameters to provide. Each c 2 
can assume any nonnegative value: c 2 = for the degenerate deter- 
ministic distribution; c 2 = k' 1 for an erlang Ek> the sum of k i.i.d. 
exponential random variables; and c 2 > 1 for mixtures of exponential 

2784 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1983 



distributions. Obviously, the difference between (2) and (1) can be 
large, so that (2) often significantly reduces the error. 

To obtain (2), we studied the GI/G/1 queue partially characterized 
by the moments of the interarrival-time and service-time distributions. 
Building on previous work by Holtzman, 30 Rolski, 31 and Eckberg, 32 we 
investigated the set of possible values of EW given the partial infor- 
mation. 33-37 When c 2 a > 1, formula (2) is always a possible value, i.e., 
there always is a GI/G/1 system with interarrival-time and service- 
time distributions having the specified moments in which (2) is correct. 
In general, (2) appears to be a reasonably typical value. 

For the single-node example just considered, the arrival process was 
a renewal process. More generally, it is natural to think of all the non- 
Poisson arrival processes in the model as renewal processes, either 
because they are initially renewal processes or because the algorithm 
can be interpreted as approximating general arrival processes by 
renewal processes. Hence, with one customer class, it is natural to 
think of the model as a generalization of the open Jackson network 
M/M/m queues to an open Jackson network of GI/G/m queues. Each 
node is approximated by a GI/G/m queue having a renewal arrival 
process independent of service times that are independent and iden- 
tically distributed with a general distribution. It is significant that 
QNA is consistent with the Jackson network theory: If there is a single 
class of customers, if all the arrival processes are Poisson, and if all 
the service-time distributions are exponential, then QNA is exact. 
However, for the general model few analytical results are available, so 
approximations are needed. 

The software package QNA has a flexible input procedure: the model 
will accept more than one kind of input (see Section II). For the 
standard input, only limited information is required. Only two param- 
eters are needed for each service -time distribution and each external 
arrival process. Also, a routing matrix is needed, which gives the 
proportion of those customers completing service at facility i that go 
next to facility j. (The algorithm is based on Markovian routing.) 
Hence, for n nodes, the input consists of n 2 + An numbers. 

There is also an alternate input by classes and routes. In this scheme 
there are different classes of customers and each class enters the 
network at a fixed node and passes through a specified sequence of 
nodes. For each class, there are two parameters characterizing its 
external arrival process and two parameters characterizing the service- 
time distribution at each node on its route. With this input by routes, 
different classes can have different service-time distributions at a 
given node and the same class can have different service- time distri- 
butions during different visits to the same node. For a class with n 
nodes on its route, the input consists of 3n + 2 numbers. (This includes 
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the n nodes on the route.) QNA analyzes this route input by aggrega- 
tion: All the classes are aggregated by QNA to convert the route input 
into the standard input. Afterwards, the special parameters of each 
class are used to describe its sojourn times. 

QNA also provides a fairly rich output. Several different congestion 
measures are calculated for each node: the traffic intensity (utiliza- 
tion), the expected number of busy servers (offered load), and the 
mean and variance of the equilibrium delay and number of customers 
present. In fact, for single-server nodes the delay distribution itself is 
described. Congestion measures for the entire network are also calcu- 
lated, under the approximation assumption that the nodes are sto- 
chastically independent given the approximate flow parameters. 
Means and variances of total service times, total delays, and total 
sojourn times (response times) are given. When the input is by routes, 
these characteristics are given for each customer class. Otherwise, 
these characteristics are given for any route requested by the user. 

A desirable feature of QNA is the structure of the calculus to 
transform the parameters to characterize the internal flows. The 
calculus is linear for each network operation, so that the parameters 
for the internal flows are determined simply by solving systems of 
linear equations. For the rates, the system of linear equations is just 
the familiar traffic rate equations occurring in the Jackson network of 
M/M/m queues. After having obtained the rates, we obtain the vari- 
ability parameters of the internal flows (the squared coefficients of 
variations) by solving another system of linear equations. As a by- 
product, the existence of a unique nonnegative solution for the flow 
parameters is trivially guaranteed. There is no guarantee that an 
iterative scheme will converge, and if it does, there is typically no 
guarantee that a solution is unique. The linearity also guarantees that 
the computation required is not great. Since there is only one linear 
equation per node in the network, QNA can be used to analyze large 
networks repeatedly at minimal cost. 

The linear calculus for transforming the variability parameters 
incorporates results of recent studies to improve the accuracy of the 
approximations. The general framework for approximating point proc- 
esses in Whitt 14 is used. Significant improvement over previous ap- 
proximation methods of this kind has been obtained by paying partic- 
ular attention to the difficult superposition operation. For superposi- 
tion, we use a modification of the hybrid procedure developed at Bell 
Laboratories by Albin. 15 - 1638 - 39 

We emphasize that QNA is approximate. In applications it is im- 
portant to validate the QNA output by comparing it with simulations 
and/or measurements. QNA is designed so that it is easy to incorporate 
improvements and it is easy to tune QNA for particular applications. 
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QNA also provides a useful framework for developing new approxi- 
mation procedures. Moreover, it is easy to use QNA in conjunction 
with other special algorithms available to analyze the nodes or the 
flows. 

The rest of this paper describes QNA in more detail. The paper is 
organized — just as the output is — according to the main steps in the 
analysis. The input is described in Section II. Section III describes the 
preliminary analysis to eliminate immediate feedback. The procedures 
to determine the internal flow parameters are contained in Section 
IV, and the procedures to calculate approximate congestion measures 
for the nodes are contained in Section V. Section VI contains the 
procedures to calculate approximate congestion measures for the net- 
work as a whole. 

In a sequel in this issue of the Journal, 40 we describe the performance 
of QNA by comparing it with simulation and other approximations of 
several networks of queues. The sequel illustrates how to apply QNA 
and demonstrates the importance of the variability parameters. 

II. THE INPUT 

In this section we describe the input options currently available for 
QNA. We anticipate more input options in the future. In Section 2.1 
we describe the standard input, which is relatively compact. In Section 
2.2 we describe a minor modification of the standard input, which 
allows for the creation or combination of customers at the nodes. For 
example, when a packet completes service at some node, it may cause 
several packets to be sent to other nodes. In Section 2.3 we describe 
an alternate input for different classes of customers having specified 
routes. We also describe the way QNA converts this input by classes 
and routes into the standard input of Section 2.1. 

2. 1 The standard input 

With the standard input, there is a single customer class and no 
creation of customers at the nodes. Any number of networks can be 
processed during a single run, so the user first specifies the number of 
networks. Then, for each network, the user specifies the number of 
nodes, and for each node the number of servers. For each node in the 
network, there are two parameters for the service-time distribution 
and two parameters for the external arrival process. Finally, there is 
a routing matrix, indicating the proportion of customers that go to 
node j from node i. Here is a list of the input data for each network 
with the notation we use: 

n = number of (internal) nodes in the network 

rrij = number of servers at node j 

Xq, = external arrival rate to node j 
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c% = variability parameter of the external arrival process to node j 

(squared coefficient of variation of the renewal interval in the 

approximating renewal process) 
tj = mean service time at node j 
cij = squared coefficient of variation of the service-time distribution 

at node j 
q t j = proportion of those customers completing service at node i 

that go next to node j. 
In matrix notation, Q ■ (q,>) is an n x n matrix and A ■ (Xo,) is an 
IX n vector. The user has the option of inputting tj or its reciprocal 
tij, the service rate at node j. (The same form must be used for all 

nodes.) 

The user need not specify the variability parameters clj and c£, in 
which case they are set equal to the default value one, corresponding 
to the M/M/l model having a Poisson arrival proccess and an expo- 
nential service-time distribution. (Again, this option applies to all 
nodes.) Alternatively, the user can specify only the service-time vari- 
ability parameters, cj, in which case either all the arrival-process 
variability parameters are automatically set equal to 1, yielding an 
M/G/l approximation for each node, or the QNA algorithm is applied. 

2.2 Creating and combining customers 

QNA has an option to allow creating or combining customers at the 
nodes following the completion of service. For example, a message 
processed at some node might cause messages to be sent to several 
other nodes. Alternatively, messages might be divided into packets 
after service at one node and then later recombined into messages 
after service at another node. In a job shop, the focus might shift back 
and forth between units and lots, e.g., at different nodes we might 
consider bottles, six-packs, cases, and even truckloads. 

With this option, the user must specify the multiplicative factor y 
of customer creation or combination at node j for each ;'. There is 
customer creation (combination) at node j if yj > 1 (yj < 1)- If 
customers are neither created nor combined, then yj = 1. If A, is the 
overall arrival rate to node j, then the departure rate, after this 
modification, is \jyj and the rate of departure from the network j is 



(i-i 



V» I 1 - 1 Qjk 

\ k=i 

When artificial nodes are used, the creation or combination can also 
be placed before service. 

To obtain our approximation formulas, we work with the following 
models of customer creation and combination. These models require 
integer values, but the approximation formulas and the QNA input do 

2788 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1983 



not. For customer creation, we replace each departure from node j 
with a batch of size yj. For customer combination, we replace yj x 
successive interdeparture intervals by a single one. From these models 
it is not difficult to calculate the impact of yj, e.g., the departure rate 
at node j is simply multiplied by yj. 

2.3 Input by classes and routes 

QNA provides the option of defining different customer classes. 
Each class has its own route or itinerary that specifies the sequence 
of nodes visited. Thus, for each class the routing is deterministic. Each 
class has an external arrival process that goes to the first node on the 
route. As usual, the external arrival process is characterized by rate 
and variability parameters. Also, each class may have its own service- 
time distribution at each node on its route. The service-time distri- 
butions can be different, not only for different classes, but also for 
different visits to the same node by the same class. These service-time 
distributions are also characterized by rate and variability parameters. 
(Alternatively, the user can elect to input the service-time parameters 
for each node. Then all classes have the same service-time distribution 
at each visit to a particular node.) 

As with the standard input, the user must specify the number of 
nodes and the number of servers at each node. Now we need the 
number of routes too. The required data are: 
n = number of nodes 
rrij = number of servers at node j 
r = number of routes. 
Here is a list of the input data for the /zth customer class of a network: 
rik = number of nodes on route k 
X* = external arrival rate of class k 
d = variability parameter of the external arrival process for class 

k 
Tikj = the ;th node visited by customer class k 
Tkj = the mean service time of class k at the jth node of its route 
Cskj = the variability parameter of the service-time distribution of 

class k at the ;'th node of its route. 
QNA converts this input by classes and routes into the standard 
input in Section 2.1. It then calculates the parameters of a typical or 
aggregate customer. Later, when computing sojourn times or response 
times of each customer class, QNA uses the service-time parameters 
of that customer class. The first version of QNA assumes as an 
approximation that each customer sees independent versions of the 
equilibrium distribution at each node. Hence, the waiting time before 
beginning service at each node is assumed to be the same for all classes 
and all visits. 
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We now indicate how QNA converts the input by classes and routes 
into the standard input of Section 2.1. For this purpose, let 1H be the 
indicator function of the set H, i.e., lH(x) = 1 if x G H and lH(x) = 
otherwise. 

First, we obtain the external arrival rates by 

r 

Ao, = I kl{k:n kl = j\, (3) 

i.e., the external arrival rate at node j, Xq,, is the sum of all route 
arrival rates X* for which the first node on the route is j. (Here the 1H 
notation is used for H = \k:n k i = ;'}.) Similarly, the flow rate from i to 
j is 

r n h -l ^ 

hj = I £ XMk, /):n k/ = i, Tiks+\ = j\ (4) 

k=l S=l 

and the flow from i out of the network is 

r 

X.o = £ X k l\k:n knk = i}. (5) 

From (4) and (5), we obtain the routing matrix Q. The proportion of 
customers that go to j from i is 



k=i 



9y= WUo+ X A* . (6) 



If node i is an active part of the network, then the denominator will 
be strictly positive. Otherwise, QNA gives an error message. 

Next, if the service-time parameters are given by routes, we obtain 
the service-time parameters for the nodes by averaging: 



Ti = 



1 £ xW{(fe, &.n k ,= j\ 

fe=l /"=! 



(7) 



The denominator in (7) will be strictly positive if node j is ever visited. 
Otherwise, as with (6), QNA supplies an error message. 

We obtain the node variability parameters c% using the property 
that the second moment of a mixture of distributions is the mixture 
of the second moments. Therefore, we have 

l 1 \kTlAclk,+ l)l\(k, /):n k ,= j\ 

rfidj +l) = *- 1/ " 1 r nk . (8) 

2 I i k l{(k,/):n k ,= j\ 

At this point, QNA has calculated enough information about the 
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standard input to compute the internal flow rates X, and the traffic 
intensities pj as described in Section 4.1, i.e., 

Pi = ^jTj/rrij. (9) 

QNA uses this information to calculate the variability parameters clj 
of the external arrival process. The hybrid approximation for super- 
position arrival processes in Section 4.2 is also used here because the 
external arrival process to node j is the superposition of the external 
arrival processes to node j from the different classes. If Xq, = 0, then 
clj does not matter and QNA sets clj = 1. Otherwise, 



clj = (1 - Wj) 



+ Wj 



S cl[\ k l\k:n kl = j] 
k=i 



£ \A\k:n kl = j] 



where 



w, 



*>j(Pj, h) = [1 + 4(1 - P j) 2 (uj - i)r\ 
Pj is the traffic intensity in (9), and 



v, = 



(10) 



(11) 



(12) 



I I X*l|fcn*i = ;) / J \AV'.n, x = ;) 

Example 1: To help fix the ideas, we consider an elementary example 
with n = 2 nodes and r = 3 routes. Let the number of servers at the 
nodes be mi = 40 and m 2 = 10. Let the route input be described by 
vectors 

(n k , Xk, cl\ n ku Tkl , cl kx ; . . . ; n kr , k , Tknk , c 2 sknh ). (13) 

Here suppose that the r vectors are: 

(2, 2, 1; 1, 1, 1; 1, 3, 3) 

(3, 3, 2; 1, 2, 0; 2, 1, 1; 1, 2, 1) 

(2, 2, 4; 2, 1, 1; 1, 2, 1). (14) 

The first route corresponds to a Poisson arrival process at rate 2 to 
node 1, with all customers being fed back immediately for a second 
service before departing from the network. (Of course, the arrival 
process need not actually be Poisson; a Poisson process always has 
c 2 = 1 but other processes could have c 2 = 1 too.) The second class 
also enters at node 1, then goes to node 2 and back to node 1 before 
departing from the network, etc. 

By (3), the external arrival rates are X i = 5 and X 02 = 2. By (4), the 
internal flow rates are X n = 2, Ai 2 = 3, X 2 i = 5, and X 22 = 0. By (5), 
the flow rates out of the network are Xi = 7 and X 20 = 0. By (6), the 
routing probabilities are: q n = 1/6, q 12 = 1/4, q 2i = 1, and q 22 = 0. By 
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(7), the mean service times are ti = 2 and r 2 = 1. By (8), c Bi = 1.67 
and c 2 2 = 1.00. Note that both service times at node 2 in (14) have 
mean 1 and squared coefficient of variation 1, as with a common 
exponential distribution, so we should want t 2 = c a2 = 1. 

To obtain the internal arrival rates, we solve the traffic rate equa- 
tions as in Section 4.1, i.e., 

n 

h = Ao, + £ KQij, ( 15 ) 

i-l 

to obtain Xi = 12, X 2 = 5, p x = 0.6, and p 2 = 0.5. 

Finally we obtain the variability parameters c<y. First, from (12), 
~ Vl = 25/13 and v 2 = 1.0. Then, from (11), ibi = 0.629 and u) 2 = 1, so 
that clx = 1.38 and c§ 2 = 4. Since there is only one external arrival 
process to node 2, we should have i/ 2 = ii) 2 = 1 and c 02 = Cs = 4. 

III. ELIMINATING IMMEDIATE FEEDBACK 

In this section we describe a function that QNA can perform before 
calculating the internal flow parameters and analyzing the congestion. 
The user can elect to reconfigure the network to eliminate immediate 
feedback. This procedure, which was originally suggested by Kuehn, 22 
usually improves the quality of approximations. Hence, it is recom- 
mended and is performed in the standard version of QNA. 

Immediate feedback occurs whenever q u > 0. Since QNA assumes 
Markovian routing, each customer completing service at node i is 
immediately fed back to node i to be served again with probability q u . 
Each time the customer goes to the end of the line. With the decom- 
position method, QNA assumes the customer finds the equilibrium 
number of customers at the node each time, with each visit being an 
independent experiment. 

QNA eliminates immediate feedback by giving each customer, upon 
arrival from another node, his or her total service time before going 
to a different node. This is equivalent to putting a customer immedi- 
ately fed back at the head of the line instead of at the end of the line. 
Transitions from node i back to node i are eliminated and the new 
probability of a transition to node j becomes the old conditional 
probability given that the customer departs from node i. In other 
words, each visit to node i from elsewhere plus all subsequent times 
immediately fed back are interpreted as a single visit. The service time 
is increased to compensate. 

The motivation for this procedure is easy to explain. For a multi- 
server node with Bernoulli (Markovian) feedback and iid service times 
that are independent of a general arrival process (not necessarily 
Poisson or renewal), the distribution of the queue length process (but 
not the waiting times) is the same after this transformation. Hence, 
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we calculate the approximate values of the mean and variance of the 
equilibrium queue length for the transformed node without feedback 
and use them to derive approximate waiting time characteristics. By 
Little's formula, 41 ' 42 the expected waiting time is also exact, i.e., the 
only error is in the approximation of the arrival process by a renewal 
process and the approximations for the characteristics of the GI/G/m 
queue; there is no additional error due to the immediate feedback. 

The first step of the reconfiguring procedure is quite simple: the 
new service time is regarded as a geometric mixture of the n-fold 
convolution of the old service-time distribution. The parameters t,, 
c'ii, and <7y are changed to f ,-, c 2 „ and q^ when qu > 0: 

f, = 17/(1 - q a ) 

cli = qa + (1 - qu)cli 
qu = 

Qij = Qs/(l ~ Qa), J * i- (16) 

Afterwards, when calculating congestion measures for node i, QNA 
makes further adjustments. When we eliminate immediate feedback 
according to (16), we no longer count the times a customer is fed back 
immediately as separate visits. Hence, we need to adjust the congestion 
measures that are expressed per visit. For example, since the expected 
number of visits to node i per visit from outside is (1 — qu)~\ to obtain 
the expected waiting time per original visit to node i, we multiply the 
values of the expected waiting time EW t obtained from (16) by 
(1 - qa). Of course, the number of customers at each node is not 
affected by the feedback treatment. 

Let \t, f ,-, etc., represent the new adjusted values. In terms of the 
parameters X,, t,, etc., obtained using (16), the new adjusted values 
are: 

X, = X,/(l - q u ) 
f, = (1 - qu)n 
cl = (cl - <7.,)/(l - qu) 
EWi = (1 - q a )EWi 
VtnW) = (1 - 9a )Var(T/) - c 2 Bl rf 
Var(T/) = c 2 (77)(W, + t,) 2 
c 2 (T!) = c 2 (f/)(l + q u ) + q u 
c 2 (T!) = (Var W( + c&fKEWi + f,)" 2 
Var(jy/) = EN.clr 2 + cHNMEN^il (17) 
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where N t represents the equilibrium number of customers at node i 
and the T, variables represent the sojourn time per visit at node i. 

We obtain the variables t, and c 2 ,- in (17) by inverting the operation 
in (16), so we receive the original data again. The last five formulas in 
(17) involving the second-moment characteristics of W, are based on 
the results of Ni in the transformed system and heavy-traffic limit 
theorems for networks of queues by Reiman. 28,29 The main quantity 
desired is Var(W,); the variable W- is a preliminary approximation 
for Wi. 

In heavy traffic, the changes in the queue length at the nodes are 
negligible during a customer's sojourn in the network. Hence, if node 
i is visited X, times by some customer, then the total sojourn timeat 
node i, say T(, is distributed approximately (in heavy traffic) as X,T/, 
where X, is independent of T[ and T[ is the sojourn time per individual 
visit in (17). (We use T{ and T[ instead of T, and T, because we do 
not use the description of T, obtained directly from (16) and T, will 
differ from T /.) By the independence, ET[ 2 = EXfET! 2 . Since X, is 
geometrically distributed with mean (1 - qu)~\ c 2 (X.) = q iit and we 
obtain the seventh formula in (17). 

The sixth and eighth formulas in (17) just express the formula for 
c 2 in terms of the mean and variance and the fact that the sojourn 
time is the sum of a waiting time and a service time. The final formula 
for Vai(W-) is obtained by approximating W{ by the sum of N t iid 
service times, using standard formulas for the variance of a random 
sum (e.g., compute EW- 2 by firs^ conditioning on N { ). Finally, we 
obtain the fifth formula for Var(tPi) by splitting the variance of 27 
into waiting-time and service-time components and dividing by the 
expected number of visits to node i. As a consequence, Var(T-) seems 
more reliable than Var( W { ). This procedure makes Var(T,), computed 
from Var(Wi) by adding variance components as in Section VI, agree 
with the direct formula for Var(T/ ) in (17). 

The congestion measures based on (16) can be used to describe the 
total delays and total sojourn times of arbitrary customers in the 
network as in Section 6.2, but the congestion measures based on (17) 
are needed to describe the behavior of particular customers with 
specified routes as in Section 6.3. However, as stated above, Var(T/ ) 
in (17) is an attractive alternative to Var(T,) obtained via (16). 

Experience indicates that eliminating immediate feedback often 
yields a better approximation (see Kuehn 22 and Sections V and VII of 
Whitt 40 ). It is also often desirable to reconfigure the network to 
eliminate almost-immediate feedback, e.g., flows that return relatively 
quickly after passing through one or more other nodes (see Section V 
of Whitt 40 ). Further study is needed to understand feedback phenom- 
ena and to develop improved approximations. 
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IV. THE INTERNAL FLOW PARAMETERS 

In this section we indicate how QNA calculates the internal flow 
parameters. In Section 4.1 we focus on the flow rates, which are 
obtained via the traffic rate equations, just as with the Markov models. 
In Section 4.2 we display the corresponding system of linear equations 
yielding the variability parameters. The remaining subsections explain 
how the variability parameter equations were obtained. The basic 
operations of superposition, splitting, and departure are discussed in 
Section 4.3, Section 4.4, and Section 4.5, and their synthesis in Section 
4.6. 

4. 1 Traffic-rate equations 

In this step QNA calculates the total arrival rate to each node. Let 
Xj be the total arrival rate to node j, let jj be the multiplicative factor 
of customer creation at node j as specified in Section 2.2, and let 5, be 
the departure rate (to other nodes as well as out of the network) at 
node j. In general, 8j = X/y/. If there is no customer creation, then 
yj = 1 and the rate in equals the rate out. 

The fundamental traffic-rate equations are just 

n 

A; = Aq, + I A,7.<7.y (18) 

i=i 

for j = 1, 2, • • • , n, or in matrix notation 

A = Ao(/ - TQ)~\ (19) 

where A — (Aq/) is the external arrival-rate vector, Q ■ (<?,;,) is the 
routing matrix, and T = (yy) is the diagonal matrix with 7,, = 7, and 
7y = for i # j. When there is no customer creation, 7, = 1 and T = 
I. Of course, (18) is just a system of linear equations. To solve them is 
equivalent to inverting the matrix (/ — TQ) in (19). When customers 
can be created at the nodes as in Section 2.2, special care should be 
taken to be sure that (18) has a solution. We need to have sp{TQ) < 1 
where sp(TQ) is the spectral radius of TQ. 

Given the arrival rates, it is possible to solve for the traffic intensities 
or utilizations at each node, defined by 

Pi = XiTi/mi, 1 « i < n. (20) 

If pi ^ 1, then the ith node is unstable. If any node is unstable, the 
algorithm gives an error message, prints out the traffic intensities, and 
stops. The associated offered load at node i, which coincides with the 
expected number of busy servers [see p. 400 of Heyman and Sobel 41 
or (4.2.3) of Franken et al. 42 ] is 

a, = X,t„ 1 ^ i =S n. (21) 
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The parameters a, and p, coincide for a single server, with a, tending 
to be more useful as the number of servers, m it increases (obviously 
when m, = o°). 

After the arrival rates have been calculated for the nodes, QNA 
calculates related quantities for the arcs: 

Xy = Xi7.<7y —the arrival rate to node j from node i 

Pu = XyA; — tne proportion of arrivals to j that 

came from i, i ^ 0. (22) 

Similarly, QNA calculates the following output rates: 

di = X,7, 11 — 2 Qu ) —the departure rate out of the 
V J =1 i network from node i 

n 

d = Y di — the total departure rate out 

1=1 of the network. (23) 

4.2 Traffic variability equations 

The heart of the approximation is the system of equations yielding 
the variability parameters for the internal flows, i.e., the squared 
coefficients of variation for the arrival processes, e%. (These are 
derived in Sections 4.3 through 4.7.) The equations are linear, of the 
form 



4 = aj+ I clibij, l^j^n, (24) 

;=i 

where a, and 6,> are constants, depending on the input data: 
aj—l + Wj j (pojclj - 1) 

+ E p y [U - in) + (i - "ij)y^ijphi)\ 



(25) 



and 



bij = WjPijqijjiivij + (1 - Pij){l - pi)), (26) 

where jc„ j/ y , and Wj depend on the basic data determined previously, 
e.g., pi, m, and c%, but not on the variability parameters 0% being 
calculated. The parameter 7, is the multiplicative factor of customer 
creation or combination, introduced in Section 2.2. The variables x t 
and vij are used to specify the departure operation; the variable Wj is 
used to specify the superposition operation. The variables vg and Wj 
are weights or probabilities that are used in convex combinations 
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arising in hybrid approximations for departure and superpositions, 
respectively. The variables xj, vtj, and Wj are included to make modifi- 
cation of the algorithm based on (24) easy. The specific values in this 
version of QNA are: 

Xi i - 1 + mr 05 (max{c 8 2 ,, 0.2} - 1), (27) 



and 



vn = 0, 






(28) 


wj = [1 + 4(1 - Pj ) 2 (uj - 


- I)]" 1 


(29) 


"j = It P§ 
L'=o . 


-i 




(30) 



with 



andpy in (22). 

It is significant that it is easy to modify this system of equations. 
For example, other hybrid procedures for departures or superpositions 
can be introduced just by changing vq and Wj, respectively. In this way, 
it is easy to calculate and compare the variability parameters for 
several different approximation procedures. 

4.3 Superposition 

The purpose of the following sections is to explain the key approx- 
imation equations (24) through (30), which yield the variability param- 
eters for the internal flows. The approximations are all based on the 
basic methods in Whitt: 14 the asymptotic method and the stationary- 
interval method. We consider the" basic operations — superposition, 
splitting, and departure — in turn, and then their synthesis. 

For superposition, the stationary-interval method is nonlinear so it 
presents difficulties. 14 " 16,22 Moreover, there appears to be no natural 
modification that makes it linear. On the other hand, the asymptotic 
method is linear. By the asymptotic method, the superposition squared 
coefficient of variation c\ as a function of component squared coeffi- 
cients of variation cf and the rates A, is just the convex combination 

d = l Uill\k\ c 2 . (31) 

However, neither the asymptotic method nor the stationary-interval 
method alone works very well over a wide range of cases, e.g., see 
Section III of Whitt. 40 Albin 1516 found that considerable improvement 
could be obtained by using a refined composite procedure, which is 
based on a convex combination of c% for the asymptotic method and 
cli for the stationary-interval method. Her hybrid ch is of the form 
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c 2 h = wcl + (1 - w)ch. (32) 

Unfortunately, since ch is nonlinear, so is c&. However, Albin found 
that a convex combination of c\ and the exponential c 2 of 1 worked 
almost as well, having 4-percent average absolute error as opposed to 
3 percent. Hence, we use such a hybrid procedure, namely, 

ch = wcl + (1 - w) 

= wl (hllXk) cj + l-w, (33) 



where w is a function of p and the rates. Extensive simulation 
prompted Albin to suggest the weighting function 

w = [1 + 2.1(1 - p)"*] -1 , (34) 

where 



v = 



k 



)T 



I A./IM . (35) 



Note that if there are k component processes with equal rates then 
v = k. The parameter v can be thought of as the number of component 
streams, with it being an equivalent number if the rates are unequal. 
However, the weighting function (58) fails to satisfy an important 
consistency condition: We should have w = 1 when v — 1; if there is a 
single arrival process, the superposition operation should leave the c 2 
parameter unchanged. Moreover, new theoretical results 39 indicate 
that the exponent of (1 - p) in (34) should be 2. Hence, we use formula 
(33) based on the weight function w in (29). 

4.4 Splitting 

No approximation is needed for splitting because a renewal process 
that is split by independent probabilities (Markovian routing) is again 
a renewal process. However, approximation is of course indirectly 
associated with this step because the real process being split is typically 
not a renewal process and the splitting is often not according to 
Markovian routing. 

Since a renewal process split according to Markovian routing is a 
renewal process, the asymptotic method and the stationary-interval 
method coincide. If a stream with a parameter c 2 is split into k streams, 
with each being selected independently according to probabilities p,, 
i = 1, 2, • • • , k, then the ith process obtained from the splitting has 
squared coefficient of variation cf given by 

c 2 = Pi c 2 + 1 - Pi , (36) 

which is clearly linear. Formula (36) is easy to obtain because the 
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renewal-interval distribution in the split stream is a geometrically 
distributed random sum of the original renewal intervals. 

4.5 Departures 

For the stationary-interval method with single-server nodes, we 
apply Marshall's formula for the squared coefficient of variation of an 
interdeparture time, say cl, in a GI/G/1 queue: 43 ' 44 

d = cl + 2p 2 cl - 2p(l - P )nEW, (37) 

where EW is the mean waiting time. Since EW appears in (37), the 
congestion at the node affects the variability of the departure process. 
A stationary-interval method approximation for cl is obtained by 
inserting an approximation for EW in a GI/G/1 queue. Our analy- 
sis 33-37 suggests that it suffices to use the linear approximation (2) 
with g set equal to one. When this is combined with (37), we obtain 
the simple formula 

d = p 2 cl + (1 - p 2 )cl (38) 

A simple extension of (38) for GI/G/m queues that is being used in 
the current version of QNA is 

d = 1 + (1 - p 2 )(cl - 1) + -S- (cl - 1). (39) 

vm 

Note that (39) agrees with (38) when m = 1 and (39) yields cl = 1 as 
it should for M/M/m and M/G/oo systems for which the stationary 
departure process is known to be Poisson. The third term in (39) 
approaches as m increases, reflecting the way multiple servers tend 
to act as a superposition operation. A basis for further refinements of 
(39) is the asymptotic analysis of departure processes in Whitt. 45 This 
asymptotic analysis shows that in some cases the variability of the 
departure process depends on the arrival and service processes in a 
more complicated way. 

As with superposition, the asymptotic method yields a more elemen- 
tary approximation than the stationary-interval method. In fact, the 
asymptotic-method approximation for the departure process is just 
the arrival process itself, i.e., the asymptotic -method approximation 
for cl is just Ca. 44 The number of departures in a long interval of time 
is just the number of arrivals minus the number in queue, and the 
number in queue fluctuates around its steady-state distribution, 
whereas the number of arrivals goes to infinity. 

It remains to combine the basic methods to form a refined hybrid 
procedure. However, limited experience indicates that this refinement 
is not as critical as for superposition. The stationary-interval method 
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alone seems to perform better for departure processes than for super- 
position processes. 44 

The most appropriate view for the departure process— the station- 
ary-interval method or the asymptotic method — depends on the traffic 
intensities at the next nodes where the departures are arrivals. As the 
traffic intensity of the next node increases, the asymptotic-method 
approximation for the departure process becomes more relevant. For 
example, consider the case of two queues in series with parameters 
Xoi, Coi, in, Csi, M2, and c» 2 . If Hz -* X while ah remains unchanged, then 
p 2 — > 1 and the second queue is in heavy traffic. Under such heavy 
traffic conditions, it has been shown 26 that the congestion measures 
at the second node are asymptotically the same as if the first facility 
were removed, i.e., as if the arrival process to the second node were 
just the arrival process to the first node. More generally, for any arrival 
process it has been proved that the asymptotic method is an asymp- 
totically correct approximation for a queue in heavy traffic. 26 

Hence, it is natural to tune the departure approximation by using 
the traffic intensities in the following nodes. Since the departure 
process typically will be split and sent to different nodes with different 
traffic intensities, it is appropriate to do the tuning after splitting. Let 
ci be the departure c 2 at node i. Then 

4 = QiA +l~Qij (40) 

is the c 2 for the portion of the departures going to node j. We let c 2 , be 
a weighted combination of the approximations obtained by the asymp- 
totic method and the stationary-interval method [using (39)]: 

4 = ViMacli + 1 - qij) 

+ (1 - "i/)[<7.y[l + (1 ~ p1)(& - 1) 

+ p 2 im- Qb (cl - 1)] + 1 - qui (41) 

where vy is chosen to satisfy ^ v% ^ 1 and be increasing in p, with 
Vjj -* 1 as pj — * 1. However, we have not yet found that positive py 
helps, 44 so the current version of the QNA uses (28). 

From (38) it is clear that the departure process variability, as 
depicted by QNA, is an appropriate weighted average of the arrival- 
process variability and the service-time variability. Hence, when the 
service time is deterministic, so that c% = 0, the departure process is 
less variable than the arrival process. However, the actual reduction 
of variability in a network caused by deterministic service times often 
is not as great as predicted by (38) or (39). Hence, we have replaced 
(39) by 
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d - 1 + (1 - P 2 )(cl - 1) + 4= (max{c 8 2 , 0.2} - 1). (42) 

vm 

After making this change, we get (25) through (28). 

4.6 Customer creation or combination 

We treat customer creation or combination as a modification of the 
departure process. When there is customer creation at node i, we 
replace each departure by a batch of size 7,. When there is combination 
at node i, we replace each interdeparture interval by the sum of 7, _1 
such intervals. These make more sense for integer values, but we do 
not require it. Hence, as described in Section 2.2, the departure rate 
from node i is 7,-X,- when the arrival rate is A,. We use the asymptotic 
method to obtain the variability parameter. Since the number of 
departures from node i in a large time interval is 7, times the number 
of arrivals, the asymptotic-method approximation of the variability 
parameter for customer creation or combination is just to multiply 
ci by 7,. (By the asymptotic method, c 2 = lim t — »Var N(t)/EN(t); see 
Section 2 of Whitt. 14 ) This is done before splitting. 

4.7 Synthesis 

We obtain the basic system of equations (24) through (30) by 
combining Sections 4.3 through 4.6 as follows: 

n 

Cli = 1 - Wj + Wj X PijC'fj 
=0 

n 

- 1 - wj + Wj pojclj + X Pij{vij[qijyiC 2 ai + (1 - qtj)] 
1=1 

+ (1 - P»){TJ0f[l + (1 - Pf)(d ~ 1) 

+ p 2 mr 05 (max|cL 0.2) - 1)] + 1 - Qij \) . (43) 

The first line is based on superposition, Section 4.3, and the second 
line is based on departure, splitting and customer creation, Sections 
4.4 through 4.6. 

V. CONGESTION AT THE NODES 

Having calculated the rate and variability parameters associated 
with each internal arrival process, we are ready to calculate the 
approximate congestion measures for each node. At this point we have 
decomposed the network into separate service facilities that are ana- 
lyzed in isolation. Each facility is a standard GI/G/m queue partially 
characterized by five parameters: the number of servers plus the first 
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two moments of the interarrival time and the first two moments of 
the service time. Instead of the moments we use the arrival rate X, the 
mean service time t, and the squared coefficients of variation c 2 and 
c 2 . Since we are focusing on a single node, we omit the subscript 
indexing the node throughout this section. 

There are many procedures that could be applied at this point. 
We could fit complete distributions to the parameters, 14 and then 
apply any existing algorithm for solving a GI/G/m queue or a special 
case. Among the attractive options are procedures for analyzing the 
GI/G/1 queue, 46 the M/PH/m queue with phase-type service-time 
distributions, 4752 the GI/H k /m queue with hyperexponential service- 
time distributions 53,54 and the GI/E k /m queue with Erlang service- 
time distributions. 65 Also available are approximations based on 
heavy-traffic and light-traffic limiting behavior. 56,57 The actual proce- 
dures used in this version of QNA, however, are quite elementary. Our 
study of the GI/G/1 queue 33-37 indicates that these elementary proce- 
dures are consistent with the limited information available. Since the 
arrival process is usually not a renewal process, and since only two 
moments are known for each distribution, there is little to be gained 
from more elaborate procedures. In fact, a user of QNA should be 
cautioned not to rely too heavily on detailed descriptions such as the 
tail of the waiting-time distribution. Such detailed descriptions may 
prove to be reasonably accurate, but they should certainly be checked 
by simulation. 

We now describe the congestion measures provided by QNA. In 
Section 5.1 we treat the single-server node and in Section 5.2 we treat 
the multiserver node. 

5.1 The Gl/C/I queue 

We begin with the steady-state waiting time (before beginning 
service), here denoted by W. The main congestion measure is the 
mean EW, but we also generate an entire probability distribution for 
W. First, the approximation formula for the mean is as in (2): 

EW = rp(cl + d)g/2(l - p), (44) 

where g = g(p, cl, cl) is defined as 



, exp 
g(p, cl, cl) = 

1. 



2(1 - p) (1 - cl) 2 
3p cl + cl . 


, C a 2 < 1 




c 2 3*l 



(45) 



When cl < 1, (44) is the Kraemer and Langenbach-Belz approxima- 
tion, 68 which is known to perform well. 33 " 37,59 When cl > 1, the original 
Kraemer and Langenbach-Belz refinement does not seem to help, so 
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it is not used. Note that (44) is exact for the M/G/l queue having 
ci-1. 

Let the number of customers in the facility, including the one in 
service, be denoted by N. The probability that the server is busy at an 
arbitrary time, P{N > 0), and the mean EN can be obtained from 
Little's formula (see Section 11.3 of Heyman and Sobel 41 ): 

P(N > 0) = p (46) 

and 

EN = p + \EW. (47) 

Formula (46) is exact even for stationary nonrenewal arrival processes 
and (47) is exact given EW. 

For the probability of delay, P(W > 0), denoted here by a, we use 
the Kraemer and Langenbach-Belz approximation: 58 

a = P( W > 0) = p + (cl - I) fi(l - p)h(p, & d), (48) 

where 

1 + cl + pel 2 

l + p(c 2 -l)+p 2 (4c 2 + c 8 2 )' Ca ^ i 

hip,d (■;;) = < (49) 

ci + pHid + ciy Ca ^ i - 

Formula (48) also yields the correct value for M/G/l systems, namely, 
p. Additional supporting evidence for (48) is contained in Whitt. 60 

We next focus on the conditional delay given that the server is busy, 
denoted by D. Obviously, ED = EW/a. We first give an approximation 
formula for the squared coefficient of variation of D, cd. This formula 
is the exact formula for the M/G/l queue, with the service-time 
distribution being H\ when c 2 5* 1 and Ek when c 2 = k~\ where H\ is 
the hyperexponential distribution with balanced means and Ek is the 
Erlang distribution (see p. 256 of Cohen 61 and Section 3 of Whitt 14 ). 
The idea underlying this approximation is that the conditional delay 
D in a GI/G/1 queue (rather than the total delay W) depends more 
on the service-time distribution than on the interarrival-time distri- 
bution. Hence, the M/G/l formula for cd is used as an approximation 
for all GI/G/1 systems. The M/G/l formula for c D is: 

cl - 2p - 1 + 4(1 - p)ds73(c 2 + I) 2 , (50) 

where d* = E(v 3 )/(Ev) 3 with v being a service-time random variable. 
Even E(v 3 ) is available, it can be used in (50), but since E(v 3 ) is not 
available with two parameters, we use approximations for dl. The 
approximations are based on the H\ and E k distributions. 
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Case 1 : When c 2 2* 1, 

dl = 3c 8 2 (l + c 2 ), (51) 

which comes from the H \ formulas: 



a * 4 



[I 1 1 

M 2 + (1 " 9) 2 J 



and 

g = [1 + V4(c 2 - l)/(c 8 2 + l)]/2. 

Case 2: When c 2 < 1, 

d 3 a = (2c 2 + l)(c s 2 + 1). (52) 

We obtain formula (52) by considering an Erlang E k variable, which 
can be represented as the sum of k iid exponential random variables 
Xi with mean r/k, where t is the mean of the E k variable. In this case 

E{X, + • • • + X k f = kE(Xl) + Sk(k - l)E{X\)E{X x ) 

+ k(k - l)(k - 2)(EX 1 ) 3 



= (^j [6k + 6k(k - 1) + k(k - l)(k - 2)] 
so that 

«-*i¥±i» = ( 1 + f)( 1+ J), 

which reduces to (52) because cl = k' 1 for an E k variable. Note that 
(51) and (52) agree at the boundary when c 2 = 1. 

From (44), (48), and (50) through (52), we immediately obtain 
formulas for Var(D) and ED 2 : 

Var(D) = (ED) 2 d = (EW)*ch/o* 

E(D 2 ) = Var(D) + (ED) 2 . (53) 

From D we then obtain second-moment characteristics for W: 

2 £(W 2 ) _ aE(D 2 ) _ cl + \-o 

Cw ~(EW) 2 (aED) 2 

Var(W) = (EW) 2 c 2 w and E(W 2 ) = Vai(W) + (EW) 2 . (54) 

We now indicate how QNA calculates an approximate probability 
distribution for W. The distribution has an atom at zero as given in 
(48) and a density above zero. The density is chosen so that W and D 
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have the first two moments already determined for them. (This is the 
general rule, but it is not quite followed in Cases 2 and 4 below.) 

Case l:ch> 1.01. Let D have the H\ density (hyperexponential with 
balanced means) 

/d(x) = P7i<?~ 711 + (1 ~ Ph*-**, x**Q, (55) 

where 



p = [1 + V(c 2 D - l)/(c 2 D + l)]/2, 
7, = 2p/ED and 72 = 2(1 - p)/ED. (56) 

Case 2: 0.99 s£ Cd *£ 1-01. Let D have the exponential density with 
mean ED. 

Case 3: 0.501 «£ Cd < 0.99. Let the distribution of D be the convo- 
lution of two exponential distributions with parameters 71 and 72 
(71 > 72), i.e., let D have density 

/b(x) = ( 7l72 ) (e^ 21 - e"W), x > 0, (57) 

\7i ~ 72/ 

where 



. £D + V2 Var(D) - (fiD) 1 

7? = 2^ 



and 



yj 1 = ED- 72 '. (58) 

The associated tail probabilities are 

P(D >x) = foe-™ - 72e- 7lX )/(7i - 72). (59) 

Case 4: c\> < 0.501. Let D have an E 2 (Erlang) distribution with 
mean ED, which has c 2 = 0.5. Its density is 

f D (x) = y 2 xe-^, x^0, (60) 

where 7 = 2/ED. The associated tail probabilities are 

P(D > x) = e~ yx (l + yx), x>0. (61) 

For deterministic service times, d 3 s = 1, so that the smallest possible 
Cd via (50) is (1 + 2p)/3. Hence, Case 4 above will not occur often. 

Finally, we come to the second moment and variance of N, the 
number in system. For the M/G/l queue, it is not difficult to compute 
E(N 2 ). Since the steady-state number in the facility is equal to the 
number of arrivals during a customer's time in the facility, it is easy 
to compute the moments of N from the moments of W; for example, 
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E(N 2 ) = X(EW + Ev) + X 2 [E( W 2 ) + 2EWEv + E(v 2 )] 

= XEW + p + X 2 E(W 2 ) + 2XpEW + p\c\ + 1), (62) 

and 

Var(iV) = XEW + p + p 2 c\ + \ 2 Vai(W). (63) 

We now modify the M/G/l formulas (62) and (63) for the GI/G/1 
queue. Let c% be defined by 

c% = Y l Y 2 /Y 2 , (64) 

where Yj is the M/G/l value of Var(JV) in (63) using (44) for EW and 
(54)andVar(W), 

y 2 = (i - p + <r)/max{(l - a + p), 0.000001} 

Y 3 = max{(p + XEW) 2 , 0.000001}, (65) 

and a is the probability of delay in (48). The maximum is used in (65) 
to avoid dividing by zero. For the M/G/l queue, Y 2 = 1; for GI/M/1 
queues, Y2 in (65) provides just the right correction, so that (64) is 
exact given the true value of a, EW and Var(W). The correction Y 2 
in (65) makes (64) too small for D/D/l queues by a factor of (1 + p)" 1 , 
but (64) is asymptotically correct in heavy traffic: Cat— >lasp— >lif 
either c\ > or c\ > 0. 
From (47) and (64) we immediately obtain 

Var(iV) = (EN) 2 c 2 N 

and 

E(N 2 ) = Var(iV) + (EN) 2 . (66) 

When there is immediate feedback at the node and it is eliminated, 
adjustments are necessary in the formulas of this section, as indicated 
in Section III. 

5.2 The GI/G/m queue 

The first congestion measures for multiserver nodes provided by 
QNA are exact. Even for nonrenewal stationary-arrival processes, the 
expected number of busy servers is just the offered load [see p. 400 of 
Heyman and Sobel 41 and (4.2.3) of Franken et al. 42 ]: 

E min\N, m\ = a = Xt (67) 

and the traffic intensity or utilization is 

p = a/m. (68) 

By Little's formula, as in (47), 

EN = a + XEW. (69) 
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QNA currently provides only a few simple approximate congestion 
measures for multiserver nodes. These are obtained by modifying the 
exact formulas for the M/M/m model. 18 Let characteristics such as 
EW(cl, cl, m) represent the characteristic as a function of the param- 
eters c 2 , c 2 , and m, and let characteristics such as EW(M/M/m) be 
the exact value for M/M/m system. A simple approximation for EW 
based on heavy-traffic limit theorems 26 ' 56,62,63 is: 

EW(cl, cl, m) = Py/ EW(M/M/m). (70) 

Formula (70) has frequently been used for M/G/m queues and is 
known to perform quite well in that case. 64-67 By virtue of heavy-traffic 
limit theorems, we know that (70) is also asymptotically correct for 
GI/G/m systems as p — > 1 for fixed m. Limited additional study 
indicates that (70) is also reasonable for moderate values of p when 
c\ 2* 0.9 and c 2 2* 0.9, or when c\ «£ 1.1 and c 2 =S 1.1. The actual 
value may be significantly smaller (larger) when c 2 < 0.9 and cl > 1.1 
(cl > 1.1 and cl < 0.9). 

The simple approximation (70) is also supplemented by simple 
approximations for the second moments of W and N. They are 
obtained from: 



and 



cfy(ci, cl, m) = cMM/M/m) 



c 2 N (cl, cl, m) = cUM/M/m). (71) 



Related second-moment characteristics are computed as in (54) and 
(66). 

More detailed and sophisticated approximations for multi- 
server nodes are being studied. As we indicated before, a variety of 
methods and algorithms can be applied given the parameters of the 
arrival process. 47 " 57 

VI. TOTAL NETWORK PERFORMANCE MEASURES 

In this section we describe the approximate congestion measures 
calculated by QNA for the network as a whole. In Section 6.1 we 
discuss congestion measures representing the system view, e.g., 
throughput and number of customers in the network; in Sections 6.2 
and 6.3 we discuss congestion measures representing the customer 
view, e.g., number of nodes visited and response times. In fact, there 
are actually two different customer views. In Section 6.2 we discuss 
the view of an arbitrary, typical, or aggregate customer; in Section 6.3 
we discuss the view of a particular customer with a specified route 
through the network. 
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6. 1 System congestion measures 

A basic total network performance measure is the throughput, which 
we define as the total external arrival rate Xo, 

Xo = \ 01 + . . . + Xo„. (72) 

When no customers are created at the nodes, the total external arrival 
rate equals the total departure rate from the network, so that there is 
little ambiguity about what we mean by throughput. However, when 
customers are created or combined at the nodes, as in Section 2.2, 
there is more than one possible interpretation. We might be interested 
in the rate at which arrivals are processed, i.e., (72). For example, the 
customers created at the nodes might be regarded only as extra work 
that must be done to serve the arrivals. On the other hand, we might 
be interested in the rate at which customers leave the network or in 
the rate of service completions. The departure rate from the network 
is 

d=idi=l Xm ( 1 - 1 Qij) (73) 

«=i «=i \ >=i / 

and the total rate of service completions is 

s = £ Si = £ A*,. (74) 

i=l i=l 

A description of the overall congestion is provided by the mean and 
variance of the number N of customers in the entire network. In 
general, 

EN= EN!+ • • • + EN n (75) 

and, as an approximation based on assuming that the nodes are 
independent, we have 

Var(iV) = VarM) + • • • + Var(NJ. (76) 

Formula (76) is valid for the Markovian models as a consequence of 
the product-form solution, but is an approximation in general. 

6.2 The experience of an aggregate customer 

When we turn to the congestion experienced by individual cus- 
tomers, there are two very different approaches. The first approach 
keeps strict adherence to the model assumptions with the standard 
input in Section 2.1, and is based on interpreting the routing matrix 
as independent probabilities (Markovian routing). This means that 
each time any customer completes service at node i, that customer 
proceeds to node j with probability qij, independent of the current state 
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and history of the network. If the network is cyclic, this means that 
every customer has positive probability of visting some nodes more 
than once. This is the perspective of an aggregate customer. It might 
be that no individual customer actually ever visits the same node more 
than once. 

If the aggregate view is desired, then the customer experience can 
be described by employing the basic theory of absorbing Markov chains 
as in Chapter III of Kemeny and Snell. 68 We can regard the external 
node as a single absorbing state to which all customers go when they 
leave the network or we can have more absorbing states, to distinguish 
between network departures from different nodes or different subsets 
of nodes. For this interpretation, the routing matrix Q is the transient 
subchain associated with the absorbing Markov chain and the inverse 
(/ — Q) -1 in (19) with T = I is the fundamental matrix of the absorbing 
chain (see p. 45 of Kemeny and Snell 68 ). Solving the traffic-rate 
equations is tantamount to solving for this fundamental matrix. 

From the fundamental matrix it is easy to calculate the moments of 
the number w# of visits to any state j starting from any state i (on an 
external arrival process). For example, En§ is just the (i, j)-th entry 
of (J — Q) _1 . It is also easy to calculate the probability of absorption 
into each of the absorbing states starting from any initial distribution. 
These various congestion measures are easily obtained working with 
n X n matrices. 68 

Suppose that we focus on an arbitrary, typical, or aggregate customer 
arriving on an external arrival process. Then that customer enters 
node i with probability Xo,/X , where X is defined in (72) and the 
expected number of visits to node i for each customer is 

EVt = h/Xo. (77) 

(We have used the fundamental matrix to get A,.) The mean of the 
time, T„ that an arbitrary customer spends in node i during his or her 
time in the network is thus 

ETi = (EVMn + EWt) (78) 

and the expected total sojourn time (time spent in the network from 
first arrival to final departure) for an arbitrary customer is thus 

ET = t ETi = £ EVi(n + EW S ). (79) 

The variance of T, is thus 

Var(Ti) = EVWaiW) + iV») + Var(Vi)(gWj + r,) 2 . (80) 

The term Var(Vf) in (80) as well as EVi is easily obtained from the 

fundamental matrix. In particular, Var(V,) = EV'j - {EV,) 2 and 
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EV 2 = I (XoAo)[F(2F dg - 1)],,, (81) 

where F is the fundamental matrix (/ - Q)~\ F dg is the n X n matrix 
with all off-diagonal entries and diagonal entries the same as F. 

To obtain an approximation for the variance of the total sojourn 
time in the network, we assume that the sojourn times at the different 
nodes are conditionally independent, given any particular routing. 
(This is not valid even for all acyclic networks of M/M/l nodes, 69 but 
is often approximately true. 7,70 ) In particular, for a customer entering 
some specified node and making Vj visits to node j, 1 < j ** n, before 
eventually leaving the network, 

T=(l % TA (82) 

V=i fc-1 / 

where T k j is the sojourn time for the fcth visit to node j. Our approxi- 
mation assumption is that the variables T h j are mutually independent 
given the vector (Vi, V 2 , • • • , V n ). 
Hence, 

(Vt > 2 
E(T 2 ) = I £ I T ki 
i=i \*=i / 

n n (V, Vi \ 

+ 2 S I EllTmi T,i) 

= l \EViE{T\i) + E[Vi(Vi - l)]E(T u ) 2 \ 

«=i 

+ 2 £ 2 E{Tu)E(T v )E{VtVi) (83) 

i=l 7=i+l 

so that 

Var(T) = l Var(T.) + 2 £ 2 £(T 11 )£(T i; )Cov(V 1 , Vj). (84) 
i=i «'=i y=i+i 

However, the current version of QNA ignores the covariance terms in 
(84) in the calculation of Var(T). 

6.3 The experience of a particular customer 

Another approach is to decouple the macroscopic and microscopic 
interpretations. This view is common in statistical mechanics. The 
total network may exhibit statistical regularity not evidenced in any 
single particle (customer). In this view, we think of the total system 
evolving as if customers were routed according to independent proba- 
bilities, even though individual customers may have very different 
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routing probabilities, perhaps nonrandom routing or acyclic routing. 
For example, we may consider the cyclic network entirely appropriate 
for the macroscopic view even though no individual customer ever 
visits any node more than once. In order for this view to be realistic, 
each individual customer should have a relatively negligible effect on 
the total network. 

The procedure here is to solve for the equilibrium or macroscopic 
behavior of the network first and then afterwards consider particular 
customers. The particular customers will have their own routes 
through the network and perhaps their own service times at the nodes 
along the way. There are two cases, depending on whether the input 
is by classes and routes, as in Section 2.3 or the standard input as in 
Section 2.1. 

6.3. 1 Input by classes and routes 

First, suppose that we are using the input by classes and routes in 
Section 2.3. Then the particular customers correspond to the customer 
classes specified in the input. Hence, each customer has a deterministic 
route through the network and possibly special service times at the 
nodes on the route. In this case, as described in Section 2.3, QNA first 
converts the input by classes and routes into the standard input in 
Section 2.1. Then QNA solves for the equilibrium behavior. Finally, 
congestion measures are calculated for the different classes under the 
assumption that they follow their originally specified special routes 
and that upon arrival at the nodes on the route they see independent 
versions of the equilibrium state of the network. Hence, in the notation 
of Section 2.3, for a customer in class k, the expected total service 
time is 

1 T kj , (85) 

y=i 

the expected total waiting time is 

I E(W nv ), (86) 

;=i 

and the expected total sojourn time or response time is the sum of 
(85) and (86). Similarly, for a customer in class k the variance of the 
total service time is 

2 rljclkj, (87) 

;'=i 

the variance of the total waiting time is 
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J Vai(W nv ), (88) 

and the variance of the total sojourn time is the sum of (87) and (88). 

6.3.2 The standard input 

With the standard input in Section 2.1, the user must specify the 
particular customers to be analyzed. In this case, the user specifies 
classes with routes and possibly service times (rate and variability 
parameters), but these data are not used in calculating the equilibrium 
behavior. The decoupling principle is used with greater force here; 
there need not be any consistency between the microscopic and mac- 
roscopic views: This additional input does not affect the equilibrium 
behavior of the total network. 

In the current version of QNA the individual customer routes are 
deterministic, so that the additional input required is just as in Section 
2.3 and the congestion measures are just as in (85) through (88) in 
Section 6.3.1. However, it is possible to modify QNA to allow random 
routes. Then the additional input would be just as in Section 2.1; for 
each class it would consist of a routing matrix plus parameters for the 
arrivals process and service times. 
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